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ABSTRACT 


Finite  difference  methods  are  used  to  solve  the  Biot  equations  for  wave 
propagation  in  a  porous  medium.  The  computational  domain  is  a  two  dimensional 
grid  of  uniform  spacing  where  truncation  of  the  grid  on  all  sides  is  accomplished  by 
applying  homogeneous  Dirichlet  boundary  conditions.  The  difference  method  is 
second  order  in  space  and  time,  and  is  seen  to  accurately  predict  phase  speeds  of  the 
primary  compressional  and  shear  waves. 
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I.  INTRODUCTION 


Mine  warfare  has  resurfaced  as  a  major  concern  for  the  United  States  after  the 
Gulf  War.  The  Tripoli  and  the  Aegis  cruiser  Princeton  suffered  mine  strikes,  driving 
home  the  danger  of  mines  and  the  need  for  effective  intelligence  and  mine 
countermeasures.  In  1992,  the  Chief  of  Naval  Operations  reorganized  the  Navy 
headquarters  establishing  a  Director,  Expeditionary  Warfare  (N85).  At  this  time,  the 
Navy  also  established  the  Program  Executive  Office  for  Mine  Warfare  -  reporting  directly 
to  the  Assistant  Secretary  of  the  Navy  for  research,  development,  and  acquisition.  [Ref.  1] 
The  United  Nations  estimates  that  during  recent  civil  and  international  strife,  more 
than  100  million  mines  have  been  laid  in  62  countries.  A  United  States  government 
report  estimates  that  the  worldwide  total  of  deployed  land  mines  increases  by  500,000  to 
1  million  each  year.  After  troops  withdraw,  land  mines  and  unexploded  ordnances 
remain  in  the  ground,  killing  and  maiming  scores  of  people  -  farmers,  travelers,  workers, 
soldiers  -  every  day  indiscriminately.  The  Secretary  general  of  the  United  Nations 
Boutros  Boutros-Ghali  now  urges  stronger  political  and  financial  support  for  efforts  to 
neutralize  mines  worldwide.  [Ref.  2] 

In  search  of  a  better  means  for  buried  mine  detection,  BBN  Systems  and 
Technology  completed  a  study  experimenting  with  Rayleigh  waves  for  mine  detection  in 
1992.  The  results  of  their  study  claimed  the  probability  of  detection  for  anti-tank  mines  is 
from  74%  to  80%  and  17  %  for  anti-personnel  mines  .  The  maximum  detection  range 
reported  was  found  to  be  15  feet.  [Ref.  3] 
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Currently,  the  Naval  Postgraduate  School  is  conducting  experiments  based  on 
Rayleigh  wave  propagation  as  mentioned  in  the  BBN  report.  The  goal  of  the  research  is 
to  improve  detection  range  and  the  probability  of  detection  of  buried  objects. 

Finite  difference  methods  applied  to  Biot  theory  is  intended  to  enhance  the 
experimental  research  mentioned.  Very  little  work  has  been  done  in  the  past  regarding 
Rayleigh  wave  propagation  in  porous  media.  Asymptotic  solutions  for  high  and  low 
frequencies  were  found  by  Deresiewicz  in  1962.  However,  the  author  has  been  unable  to 
find  numerical  experiments  on  Rayleigh  wave  propagation  in  porous  medium  valid  for 
intermediate  frequency  ranges. 

The  original  goal  of  this  thesis  was  to  produce  a  finite  difference  code  capable  of 
modeling  wave  propagation  in  porous  media  excited  by  surface  tractions.  Due  to 
numerical  instabilities  in  discretization  of  the  free  surface  boundary  conditions,  the 
problem  modeled  was  amended  to  consider  body  wave  propagation  in  porous  media.  A 
two  dimensional  solution  to  the  problem  is  discretized  using  an  evenly  spaced  grid  for  the 
"X"  and  "Y"  coordinates,  whereby  the  grid  is  truncated  by  an  application  of  homogeneous 
Dirichlet  conditions  "relatively"  far  from  the  point  of  application  of  an  incident 
displacement.  Because  the  spatial  discretization  is  uniform,  and  center  differencing  is 
used  throughout,  the  numerical  scheme  is  second  order  in  space  and  time.  The  code  is 
written  in  C++. 
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II.  THEORY 


Biot  theory  is  one  of  two  popular  theoretical  methods  to  model  wave  propagation 
in  porous  media,  the  other  theory  assumes  the  media  is  viscoelastic.  Unlike  viscoelastic 
theory,  which  is  based  more  upon  elasticity  theory,  Biot  theory  predicts  the  existence  of 
three  kinds  of  body  waves,  two  dilatational  and  one  rotational.  The  first  type  of 
dilatational  wave  involves  motions  of  the  skeletal  frame  and  the  interstitial  fluid  which 
are  nearly  in  phase  and  for  which  the  attenuation  owing  to  viscous  losses  is  relatively 
small.  The  second  type  of  dilatational  wave  is  such  that  the  frame  and  fluid  components 
are  moving  largely  out  of  phase.  Biot  derived  two  coupled  linear  differential  equations 
governing  the  relationship  between  the  displacements  of  the  fluid  and  the  frame.  [Ref.  4] 

A.  WAVE  EQUATIONS 

The  Biot  equations  are  [Ref.  4]: 

P  "  0  +  P  .i|r  =  (P  ■ -  (O' V(V  •  u) + nV2u  +  QV(V  .  U) + bF(o»(f  - f )  ( 1) 

p  +  P22^  =  QV(V  .  u)+  RV(V  •  U)  -bF(«»<f  - f )  (2) 

where 

u  =  the  displacement  vector  of  the  solid  material; 

U  =  the  displacement  vector  of  the  fluid; 

P  =  the  ratio  of  the  volume  of  the  pores  to  the  total  volume  of  the  media; 
p  =  the  shear  modulus  of  the  solid; 

Kb  =  the  bulk  modulus  of  the  free-draining  frame; 

Kf  =  the  bulk  modulus  of  the  fluid; 
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Kr  =  the  bulk  modulus  of  the  solid; 
pf  =  the  mass  density  of  the  pore  fluid; 

pr  =  the  density  of  solid  material; 
r|  =  fluid  viscosity; 

K  =  permeability; 
a  =  added  mass  term; 

F(to)=  frequency  dependent  function  to  account  for  high  frequency  attenuation; 
such  that, 

Q  =  [p(l  —  P)K,  -  PKb]  +  [1  -  P  +  pg  -  £] 

R  =  [p2Kr]*[l-p  +  p£-§] 

b  =  [T|P2]  H-  [K] 

P=[(l-P)((l-P)-§)Kr  +  p^]^[l-P  +  Pf;-^]  +  fn 

p  12  =  PfP(l  -  a) 

Pll  =  (1  — P)pr  — Pl2 
P22  =  PfP  P 1 2 

Mathematical  notation:  The  bold  letters  denote  vectors  and  the  subscripts  denote 
derivatives  with  respect  to  the  superscripted  variable. 


U= 


\ 

J 


11  = 


J 


V(V«U)  = 


"Uxx  +  Vyx  ^ 
UXy  +  Vyy  J 


V(V  •  u) 


(  ,  \ 

UxX  VyX 

V  Uxy  +  Vyy  J 


V2u  = 


UxX  +  Uyy 

^  VXX  +  Vyy 


\ 

/ 
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The  u  and  U  terms  in  equations  (.1)  and  (2)  can  be  isolated  as  follows: 

Utt(p22  -  {&)  =  QV(V  •  u)  +  RV(V  •  U)  -  bF(a)(Ut  -  ut)  -  { 

[(P  -  |i)V(V  •  u)  +  |iV2u  +  QV(V  •  U)  +  bF(co)(U,  -  ut)]};  (3) 

u,t(pi2  -  *f^)  =  QV(V  .  u)  +  RV(V  •  U)  -  bF(co)(U,  -  ut)  -  { g 

[(P  -  |i)V(V  •  u)  +  |iV2u  +  QV(V  •  U)  +  bF(co)(U,  -  u,)]};  (4) 

Rewriting  equations  (3)  and  (4)  in  scalar  form: 

Utt^P22  —  £77  =  Q(uxx  +  vyx)  +  R(UXX  +  Vyx)  -  bF(a>)(Ut  -  ut)  -  {^7 

[(P  -  (I)(UXX  +  Vyx)  +  |X(Uxx  +  Uyy)  +  Q(UXX  +  Vyx)  +  bF(tO)(Ut  -  Ut)]}  (5) 

Vu^p22  —  £77^)  =  Q(u  xy  +  Vyy)  +  R(Uxy  +  Vyy)  —  bF(co)(V  t  -  vt)  -  {£77 

[(P  -  p,)(Uxy  +  Vyy  )  +  fi.(Vxx  +  Vyy  )  +  Q(U Xy  +  Vyy  )  +  bF((0)(V t  “  Vt)]  }  (6) 
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Utt(P  12  -  Pp1P211)  =  Q(UXX  +  Vyx)  +  R(UXX  +  Vyx)-bF(G))(Ut  ut)  {Pl2 


[(P  -  (J.)(uxx  +  Vyx)  +  |i(UXX  +  Uyy)  +  Q(UXX  +  VyX)  +  bF(tO)(Ut  ~  Ut)]}  (7) 

V,,(PT2  -  f-f, r)  =  QKy  +  Vyy)  +  R(U,»  +  V„)  -  bF(»)(V,  -  V,)  ~  {g 

[(P  -  H)(Uxy  +  Vyy )  +  P(VXx  +  Vyy)  +  Q(Uxy  +  Vyy )  +  bF(tO)(V t  -  V,)]  }  (8) 

B.  BOUNDARY  CONDITIONS 

A  rectangular  grid  is  used  in  the  X  and  Y  direction.  Displacements  u  and  U  are 
set  to  zero  at  all  sides  of  the  numerical  domain.  This  creates  reflections  as  body  waves 
strike  any  of  the  sides.  Originally,  on  one  of  the  surfaces,  a  traction  free  condition  was 
applied  which  is  necessary  for  the  propagation  of  surface  waves  and  for  proper  modeling 
of  a  porous  half  space.  On  this  surface,  both  the  stresses  and  the  fluid  pressure  vanish, 
resulting  in  the  following  boundary  conditions: 


Gyyly=0  =  (P  “  2N)(UX)  +  PVy  +  Q(UX  +  Vy)l  y=0  =  f(X,  t);  (9) 

Cxyly=o  =  N(Uy+vx)ly=o  =  0;  (10) 

sly=0  =  Q(ux+Vy)  +  R(Ux  +  Vy)ly=o  =  -pxf(x,t);  (11) 

Uy  +  Vx  =  0;  (12) 
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where  o  is  the  stress  tensor  and  S  is  related  to  the  fluid  pressure  ( — (3P) ,  and  f(x,t)  is  an 
applied  normal  force  on  the  surface  of  the  half  space.  Equations  (10)  and  (12)  are  the 

conditions  that  on  the  free  surface,  shear  stress  is  zero  for  both  the  solid  and  the  fluid. 

C.  SOURCE  FUNCTION 

The  source  function  used  for  the  half  space  problem  was: 

F(x,  t)  =  Ae_(x_Xo)2e_(t_to)2 ;  (13) 

The  source  function  is  Gaussian  in  shape  along  the  x  grid  on  the  surface.  x0  is  the 
midpoint  of  the  x  grid,  while  t0  is  the  start-up  time  for  the  numerical  calculations.  'A'  is  a 
constant  set  equal  to  one  pascal. 

When  it  was  determined  that  the  discretization  for  the  free  surface  conditions  were 
unstable,  the  following  two  initial  states  of  the  displacements  were  employed  to  check  the 
accuracy  and  stability  of  the  code  at  interior  points  of  the  domain. 


i-i°  p,-Y[(>-'o)2+G-jo)2] 

J(i-io)2-Hj-jo)2 

(14) 

■>~j°  fi-Y((i-i0)2+(Ho)2] 

V(i i->o)2+(H„)2 

(15) 

.  j~j“ _  f>-Y[('-'o)2+(Ho)2] 

7<i-io)2+0-jo)2 

(16) 

i-i»  -  c-7[(i-in)2+(i-i<>)2l 

7(i-i0)2+(j-jo)2 

(17) 
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Equation  (14)  and  (15)  were  used  to  check  the  compressional  waves.  They  can  be 
seen  as  an  initial  expansion  of  the  solid  portion  of  the  porous  medium  with  two 
dimensional  gaussian  shape.  Equation  (16)  and  (17)  are  the  cross  product  of  equations 
(14)  and  (15)  with  an  out  of  the  page  vector  (k).  They  represent  an  infinitesimal  rotation 
of  the  solid  elements  of  the  two  dimensional  gaussian  distribution  of  magnitudes,  and  are 
used  to  check  the  shear  wave  of  the  porous  medium.  Both  i0  and  j0  are  set  at  the  center 
of  the  grid,  and  gamma  is  set  equal  to  0.1. 
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III.  PROGRAM  ALGORITHM 


The  program  is  written  in  C++.  The  spatial  grid  for  both  the  X  and  the  Y  axis  are 
chosen  to  be  0. 1  meter,  while  the  time  interval,  At  is  0. 1  second.  The  time  interval  meets 
the  stability  requirements  to  be  stated  shortly.  The  Biot  equations  are  scaled  prior  to  the 

running  of  the  program.  Three  different  sets  of  Biot  parameters  were  used.  Table  1 

shows  the  inputs  which  remain  unchanged  for  all  three  parameters  sets  used. 


(  . . . .  \ 

Omega  (fo)  2000  rad/s  1 

F(w) 

1  1 

Grid  spacing  (h) 

01  1 

Time  interval  (delta  ti 

_QJ_  1 

Table  1.  Input  Parameters 

1 

A.  FINITE  DIFFERENCE  SCHEME 

In  all  the  finite  difference  approximations,  center  differences  are  used. 
Nodal  points  are  found  at  (  Xi,yj  )  where: 

X;  =  iAx;  (the  horizontal  component) 

yj  =  jAy;  (the  vertical  component) 

and  the  discrete  times  at  which  the  displacements  are  found  are  tn  =  nAt; 
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Some  of  the  difference  formula's  used  are: 

uxy  =  4^l[Mi+W+l  ~  Mi+U-1  ~  +  )’ 

The  "O"  terms  are  neglected  leading  to  a  second  order  truncation  errors  throughout  the 
code. 
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When  equations  (5)-(8)  are  differenced  in  time  and  space,  we  may  then  isolate  the 


unknown  displacements  at  time  level  'n+1'  as  follows: 
Ally-1  +  Cuy-1  =  E  (18) 

BUy 1  -Duy 1  =  F  (19) 

AVy+1  +  CVy+1  =  G  (20) 

BVy1  -  Dv"|*  =  H  (21) 


The  goal  is  to  solve  for  each  grid  point's  value  at  the  next  time  level  (n+1),  to  update  (or 
shift  time  levels),  and  march  forward  in  steps  of  At. 

Solving  for  the  'n+1'  terms  from  equations  (18)  through  (21)  we  have: 


t  rn+1  _  DE+CF 

UiJ  ~  AD+BC 


(22) 


n+1  BE-AF 
Uij  AD+BC 


(23) 


n+ 1  _  DG+CH 

V  ij  "  AD+BC 


(24) 


n+1  BG-AH 
V>.j  AD+BC 


(25) 
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The  variables  A,  B,  C,  D,  E,  F,  G,  H  in  their  respective  finite  difference  form  are  as 


follows: 


A  =  bF(ra)^  +  2fbF(0))5i;;  (26) 

B  =  JT(pa-f|)+bF«»)^  +  SbF«D)5i;;  (27) 

C  =  i(p,2-£jfi)-bF(fi»i-gbF(0»5i;;  (28) 

D  =  bF(ffl)^  +  j£bF(<o)i;  (29) 

E  =  Q{p(C,J-2u!'J +<_,.,)+ 


i(vwj*i  -v"-ijn+vwj-i -vw.nj}+ 
R{^(u.>ij  -2U"j  +  V-ij  )+ 

i(v!Vlj,l-v”_,0+l+vr_lj_,-v|Vl.j_,)}- 

bF«‘»±W;'-u!7,>- 


i(Vi+l,j+l  -VP-1J+,  +v"-lj-l -Vi+1j-l  )}- 

&{£(u"+l  j  -  4<j + <-lj + <j+l + u"j-l )  }- 
-2^  +  11"^  )+ 

i(vr+1,j+.  -vil1j+i +vr_1,H  -  v^i)}- 

gbFCto^u^-U^}- 

^(pi2-^)(u^-2u^)  (30) 

F=Q{^(ur+1>j-2<j+ur_1>j)+ 

i(vf+l,j+l  ^S-lj+I  "'’VjLjj-,  -  Vi+i,j_i)}+ 

R{^(Ul+i  ,j  ~  2U"j  +  U|Lj  j)+ 

i(vr+1J+1  -  V"_j  j+1  +  VP-U-, +vp+i,j_i)}- 
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^(P-nKi^u-H  +  uf.,^ 

-  vj.,^1  +v:_,j_,  -vnWj_,  )}- 

&{^(uWo-4u"j+u"-lj  +  u!'j*l+u“H)>- 

fifQ{p(ufn  j  -  2U,"j + uf-.j  )+ 

s(p^-?l)(u!71-2u!'j);  <31> 

g = Q{i(<lti+,  -  ur.1jtl +ur_,j.,  -  u“+IJ_i )+ 

^(vy*, -2v?j+v:j_1)}+ 

R^u"^,  +ut,j_1-u|V1J_l)+ 

^(v^-av^+v^,)}- 

ww^j^r'-vr;1)- 
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K(P-t*Hi(uMo+l -“>-1.1+1  -“"+1.1-1  j+ 

i?(vU+i  -2v?j  +  vfH  )}- 

j  -  4v^  +  v"_,j  +  v"^,  +  vPj-,  )}- 

?fQ{i(uw.H  -u”-.j+. +u”-ij-i  -u:+i.h)+ 

?(v^,-2v;j+v:j_,)}- 

gbF(a.)(^)(vH-V?-‘)- 

i(p.^-£S!1)(<I,-2vrj);  (32) 

H  =  Q{i(ul,  j+,  -  u".M+l  +  u"_,  j_,  -  uf+,  j_,  )+ 
?(v!!j„-2vj+vftj_1j}+ 

R{i(uh,j+,  -u"-ij+i +ur-l<H-Uwj-,)+ 
j(v"j+,-2Vfij  +  V"j_l))- 
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g(p-w(i  k+i,j+1  -<-10+1 +<-io-i  -<+u-i  J+ 


i(v"j+I-2v1t;j+v|:H)}- 

&{£(v*ij  -4v"j  +vf_jj  +v"j+1  +  vPj.,  )}- 

pitQ{i(u.n+u+»  -ui-ijn+ui-ij-i  -u"+ij-, )+ 

i^^ij+i  —  2V|j  +  V"j_,  j}- 


Second  order  finite  difference  forms  of  the  boundary  conditions  equations  (9) 


through  (12)  are  given  by: 


(P-2N){i«+1J 


<u)}+p{i 


Q{i(ur+]j  -  UlL,  J + -  V^, )}  =  ff ; 


(34) 
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(35) 


N{i(urotl-u"j_l+v?4lJ-v;_lj)}  =  0; 

Q!i(urtij  -  of-io + v;J4,  -  vj;., )+ 

R{s(ur+,j  -  +  vrj+,  -  V”h  )  =  -pf;-,  (36) 


2 


+  Vn  ,  • 

T  v  1+1. J 


=  0; 


(37) 


Equations  (34)  through  (37)  give  the  means  to  calculate  the  'j-1'  terms  exterior  to 
the  boundary  surface  (see  Figure  (2)),  which  are  in  turn  substituted  into  equations  (30 
through  (33)  to  find  the  values  of  u  and  U  along  the  free  surface  boundary. 


B.  SCALING 

Only  two  parameters  need  to  be  scaled  in  the  Biot  equations,  time  and  length.  The 
equations  for  the  displacements  are  linear,  therefore  their  scaling  is  determined  by  the 
magnitude  of  the  applied  force  on  the  free  surface  or  by  the  initial  displacements 
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assigned.  The  time  scale  chosen  is  omega  in  rad/sec,  while  the  length  scale  used  is  L=1 
meter.  These  scales  lead  to  the  following  relationships: 

Xs  =  l;  x  =  Lxs;  <k  =  L<kI; 


ts  =  (Dt; 


t  =  — ■ 

L  GO  ’ 


d  _  „  a 
at  “  “at,  • 


(39) 


In  addition  to  scaling,  a  constant  term  of  y  is  multiplied  across  each  equation. 
The  result  is  a  constant  multiplier  ^y-  to  all  terms  involving  second  derivatives  with 
respect  to  time,  y-  to  all  first  derivatives  in  time,  and  jj  to  all  second  derivatives  in  x  or 

y- 

C.  STABILITY 

The  system  of  equations  (18)  to  (21)  was  found  to  be  stable  provided  that 


At  <  h  -s-  pc  +  Vj  ^ 2 ; 


(40) 


(41) 


(42) 


Kc  =  bulk  modulus  of  the  saturated  medium 
=  PKf  +  (l-P)Kr; 
p  =  |3pf+(l-p)pr; 

p=  shear  modulus  of  the  frame; 

h=  the  grid  interval  in  both  x  and  the  y  direction; 

At=  the  time  step; 
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VPc  =the  saturated  velocity  of  the  fast  compressional  wave; 

Vs  =the  shear  wave  velocity; 

Equation  (40)  is  the  stability  condition  used  in  finite-difference  algorithms 
involving  wave  propagation  for  elastic  media.  [Ref.  5]  For  porous  medium,  equation 
(40)  is  also  valid  because  the  saturated  velocity  of  the  fast  compressional  velocity  in  the 
porous  elastic  solid  is  normally  much  less  than  the  compressional  velocity  of  the  purely 
elastic  solid. 

After  scaling,  the  stability  equation  (40)  becomes: 

MS<(M)+  J  VI  +  V]  ■  (43) 

D.  PARAMETERS 

The  parameters  used  are  shown  in  the  table  below. 


Units 

Stoll  &  Kan 

ARL:UT  , 

ARL:UT  mod 

Bulk  parameters 

porosity 

0.47 

0.4 

0.2 

grain  density 

kg/ 171*3 

2650 

2650 

2400 

liquid  density 

kg/m*3 

1000 

1000 

1000 

grain  bulk  modulus 

Pa 

3.6E+10 

7E+09 

7E+09 

liquid  bulk  modulus 

Pa 

2E+09 

2.25E+09 

2.25E+09 

Liquid  motion 

viscosity 

kg/m  s 

IE-03 

IE-03 

0 

permeability 

m*2 

IE-10 

4.99E-11 

4.99E-11 

pore  size 

m 

IE-05 

4.99E-05 

4.99E-05 

virtual  mass  constant 

1.25 

1.75 

1.75 

Frame  response 

frame  shear  modulus 

Pa 

2.61  E+07 

2.61  E+07 

2.61  E+07 

shear  log  decrement 

... 

0.15 

0.15 

0.15 

frame  bulk  modulus 

Pa 

4.36E+07 

5.3E+09 

5.3E+09 

bulk  loa  decrement 

0.15 

0.15 

0.15 

Table  2.  Biot  model  parameters  "From  Ref.  5  and  6“ 
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The  Stoll  and  Kan  parameters  and  the  ARL:UT  parameters  are  taken  from  [Ref.  6], 
ARL:UT  parameters  are  believed  to  be  the  better  parameters  to  use  based  on  the 
experiments  performed  by  Chotiros,  ARL  University  of  Texas  at  Austin.  The  modified 
ARL:UT  is  the  same  as  ARL:UT  parameters  except  the  fluid  viscosity  is  set  to  zero, 
eliminating  any  damping  in  the  system,  and  smaller  grain  density  and  lower  porosity 
ratios  are  used.  The  wave  speeds  calculated  using  the  above  parameters  are  shown  in 
Table  3. 


f . . 

X 

Stoll  &  Kan 

ARL:UT 

ARLUT  mod 

Fast  P  wave  speed 
Slow  P  waw  speed 
Sheer  speed 

1485.14m/s 
56.21  m/s 
119.24m/s 

1694.8m/s 
192.83  m/s 
114.7  m/s 

1676.7  m/s 

1066.8  m/s 
114.07  m/s 

Table  3.  Wave  Speeds 

Equations  used  to  calculate  the  compressional  and  shear  wave  speeds  are  [Ref.  7]: 
Pfast  =  VReal((2AH)/(B  -jf-Del))  ; 


Psiow  =  7Real((2AH)/(B  -  if  +  Del))  ; 


S=  7Real(^((p22/p)-jf)/(p(C-jf)))  ; 


where, 

H  =  P  +  2Q  +  R; 

A  =  (PR-Q2)/H2  ; 

B  =  (pnR  +  p22P-2pi2Q)/(pH); 

C  =  (pnp22  — p22)/p2; 

Del=  V(B— jf)2-4A(C-jf); 
f  =  p(3/(Kpa)). 
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E.  DISCUSSION 


The  source  was  applied  to  the  finite  difference  model  only  once  at  the  initial  time 
level  in  an  attempt  to  model  an  instantaneous  point  force  surface  traction  on  a  porous 
half-space.  Because  of  the  porous  parameters  used,  the  shear  speed  was  small  compared 
to  the  fast  compressional  wave.  The  large  difference  in  speeds  is  assumed  to  be  the  cause 
of  a  numerical  instability  at  the  free  surface  boundary.  An  instability  of  this  form  was 
reported  in  the  work  of  Ilan  and  Loewenthal  in  1976  [Ref.  8].  To  test  this  hypothesis  the 
value  of  the  shear  modulus  was  increased  to  within  the  reported  stability  region  of 
center-differencing  of  the  free  surface  conditions,  but  still  the  code  was  unstable. 
Alternative  differencing  of  the  free  surface  conditions  proved  futile,  so  a  decision  was 
made  to  drop  the  free  surface  altogether  and  apply  homogeneous  Dirichlet  Boundary 
Condition  on  all  sides  of  the  numerical  domain.  This  numerical  problem  identifies  an 
area  of  future  research. 

The  Stoll  &  Kan  parameters  demonstrate  short  propagation  distances  due  to  the 
relatively  large  frequency  and  high  attenuation  of  the  media.  (Figure  (3))  The  ARL:UT 
which  neglects  fluid  viscosity  gives  a  much  better  contrast  for  the  fast  compressional 
wave  propagation  as  reflections  off  the  sides  can  also  be  observed.  (Figures  (4)  and  (5)) 
The  fast  and  the  slow  compressional  wave  are  visible  using  the  ARL:UT  mod  parameters. 
(Figures  (6)  and  (7))  From  Figure  7,  the  fast  and  slow  compressional  wave  speeds  can  be 
calculated  using  equation  (44). 
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time  steps  ^me  steps 


Solid  displacement,  ARLUT 


surface  in  x  direction 


Fluid  displacement  ARLiUT 


10  20  dO  4U  ou  bU 

surface  in  x  direction 


Figure  5.  Top  View  -  Solid  and  Fluid  Displacement,  ARL:UT 
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200 


Solid  displacement, ARL:IIT  mod 


10  20  30  40  50  60  70  80  90  100 

surface  in  x  direction 

Fluid  displacement, ARLUT  mod 


10  20  30  40  50  60  70  80  90  100 

surface  in  x  direction 

Figure  7.  Top  View  -  Solid  and  Fluid  Displacement,  ARL:UT  Mod 
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to  =  2000  rad/sec; 


L=1  m. 


(44) 


C  = 


oA  h  . 
A?  ’ 


For  the  fast  compressional  wave  speed,  it  traveled  50  grids  in  about  60  time  level  giving 
us  a  speed  of  1666  m/s.  The  slow  compressional  wave  speed  traveled  50  grids  in  100  time 
level  giving  us  a  speed  of  1000  m/s. 

To  check  the  accuracy  and  stability  of  the  code  at  interior  points  of  the  domain, 
equations  (14)  through  (17)  were  used  and  the  results  shown  in  Figure  (8)  and  (9).  In 
Figure  (8),  the  cylindrical  spreading  show  a  propagation  of  approximately  80  grids  in  100 
time  level  for  a  1600  m/s  fast  compressional  wave  speed.  In  Figure  (9),  no  shear 
propagation  in  the  fluid  is  seen  as  expected  for  fluids  do  not  support  shear  propagation. 
For  the  solid,  shear  speed  propagated  approximately  1 1  grids  in  200  time  level  giving  us  a 
shear  wave  speed  of  1 10  m/s.  It  is  extremely  difficult  to  distinguish  the  slow  wave 
propagation  from  that  of  the  fast  wave  in  Figure  (8).  However,  due  to  the  fact  that  the 
frame  and  fluid  components  of  the  slow  wave  are  moving  largely  out  of  phase  and  the  fast 
wave  moves  nearly  in  phase,  a  plot  of  relative  speed  between  the  frame  and  the  fluid 
would  make  the  slow  wave  more  visible.  (Figure  (10)).  In  Figure  (10),  the  slow  wave  is 
seen  to  propagate  to  about  50  grids  in  100  time  levels  for  a  slow  speed  of  1000  m/s. 

All  of  the  figures  shown  and  numbers  calculated  have  propagation  speeds  that 
match  well  with  the  theory  shown  in  Table  (3). 
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y  grid 


Solid, time  level  200,  source  at  center,  ARL:UT  Mod 


Fluid,  time  level  200,  source  at  center,  ARL:UT  Mod 


Figure  9.  Shear  Solid  and  Fluid  Displacement,  ARL:UT  Mod 
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IV.  CONCLUSION 


The  propagation  of  the  fast  compressional  waves  can  be  clearly  seen  in  all  three 
parameter  sets  used.  The  program  behaved  as  expected  concerning  the  body  waves  and 
their  speeds  matched  well  with  theory.  The  different  parameters  used  have  shown  the 
sensitivity  of  Biot  parameters  in  wave  propagation. 

One  area  that  should  be  pursued  is  to  determine  a  stable  free  surface  boundary 
condition  difference  operator  so  that  Rayleigh  wave  propagation  may  be  modeled. 
Another  improvement  on  the  current  code  would  be  to  include  absorbing  boundary 
conditions  on  the  sides  and  the  bottom,  to  better  model  a  semi-infinite  domain. 
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APPENDIX  PROGRAM  CODE 


y* **************************************************/ 


/*  Thesis  Program  for  Wave  Propagation  */ 

/*  Programmer:  Jonah  Shen  */ 

/*  Date:  June  1995  */ 

/*  Language:  C,  beta,  scaled  */ 


j* **************************************** **********/ 

#include  <stdio.h> 

#include  <math.h> 

/*  Initialization  of  variables  and  constants  */ 


float  beta  =  0.47; 

float  Mu  =  2.6  le7; 
float  rhor  =  2650.0; 
float  rhof  =  1000.0; 
float  Kb  =  4.36e7; 
float  Kr  =  3.6el0; 
float  Kf  =  2.0e9; 
float  N  =  le-3; 
float  K  =  le-10; 
float  delT  =  0.1; 
float  Fw  =  1.0; 
float  h  =  0. 1 ; 
float  alpha  =3.0; 
float  w  =  2000.0; 

float  a,Q,R,P,b,rhol2,rholl,rho22,AA,BB,CC,DD,EE,FF,GG,HH,II,JJ,KK,LL; 
float  MM,NN,y,f; 

float  Aa,Bb,Cc,Dd,Ee,Ff,Gg,Hh,Ii,Jj,Kk,Ll,Mm,Nn,Oo,Pp,Qq,Rr,ans; 
float  aA,bB,cC,dD,eE,fF,gG,hH,iI,jJ,kK,lL,mM,nN,oO,pP,qQ,rR; 
int  n,i,j,ind,time,tt; 


/*  beta  is  the  ratio  of  the  vol  of  */ 

/*  the  pore  to  total  vol  of  element  */ 
/*  shear  modulaus.  Pa  */ 

/*  grain  density,  Kg/mA3  */ 

/*  fluid  density,  Kg/mA3  */ 

/*  bulk  modulus  of  frame.  Pa  */ 

/*  bulk  modulus  of  solid.  Pa  */ 

/*  bulk  modulus  of  fluid.  Pa  */ 

/*  liquid  viscosity,  kg/m  s  */ 

/*  liquid  permeablility,  mA2  */ 

/*  delta  t  */ 

/*  forcing  function  */ 

/*  step  size  */ 


/*  Defining  Vectors  for  finite  difference  */ 

/*  Lower  case  is  displacement  of  points  in  the  skeletal  frame  */ 

/*  Upper  case  is  dispalcement  of  fluid  */ 

/*  Each  letter  will  be  follow  by  3  numbers,  the  first  number  */ 

/*  represent  time  second  represent  x  direction,  third  represent  */ 


33 


/*  y  direction  *! 


float  u[3][51][51]; 
float  v[3][51][51]; 
float  U[3][51][51]; 
float  V[3][51][51]; 
float  sol[300][51]; 
float  flu[300][51]; 

/*  Declare  Function  Prototype  */ 

float  FEE(float  Aa, float  Bb,float  Cc, float  Dd,  float  Ee,float  Ff, float  Gg, 
float  Hh,float  Ii,float  Jj, float  Kk,float  LI, float  Mm, float  Nn, 
float  Oo, float  Pp, float  Qq,  float  Rr); 
float  FFF(float  Aa, float  Bb,float  Cc, float  Dd, float  Ee,float  Ff,float  Gg, 
float  Hh, float  Ii,float  Jj, float  Kk,float  LI, float  Mm,float  Nn, 
float  Oo,float  Pp, float  Qq, float  Rr); 

float  FGG(float  aA,float  bB,float  cC, float  dD,float  eE,float  fF,float  gG, 
float  hH,float  iI,floatjJ, float  kK,float  lL,float  mM,float  nN, 
float  oO, float  pP,float  qQ,float  rR); 

float  FHH(float  aA, float  bB, float  cC, float  dD, float  eE, float  fF, float  gG, 
float  hH, float  il, float  jJ,float  kK,float  1L, float  mM,float  nN, 
float  oO,  float  pP, float  qQ,float  rR); 


main() 

{ 

a  =  ( 1  -beta+(beta*Kr/Kf)-(Kb/Kr)); 

Q  =  (beta*(l-beta)*Kr-beta*Kb)/a; 

R  =  beta*beta*Kr/a; 

P  =  (( 1  -beta)*((  1  -beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a  +  4*Mu/3; 

b  =  N*beta*beta/K; 

rhol2  =  rhof*beta*(  1  -alpha); 

rholl  =  ( 1  -beta)  *  rhor-rho  1 2 ; 

rho22  =  rhof*beta-rhol2; 

AA  =  b*Fw*w/(2*delT*Mu)  +  rho22*b*Fw*w/(rhol2*2*delT*Mu); 

BB  =  (w*  w/(delT*delT*Mu))*(rho22-(rho  1 2*rho  1 2/rho  1  l))+b*Fw*w/(2*delT*Mu)+ 
b*Fw*w*rhol2/(2*delT*Mu*rhol  1); 

CC  =  (w*w/(delT*delT*Mu))*(rhol2-(rho22*rhol  l/rhol2))-b*Fw*w/(2*delT*Mu)- 
b*Fw*w*rho22/(2*delT*Mu*rhol2); 

DD  =  b*Fw*w/(2*delT*Mu)+rhol2*b*Fw*w/(rhol  l*2*delT*Mu); 

/*  Initialization  of  the  arrays  */ 
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for  (n=0;n<=2;n++) 

{  for  (i=0;i<=50;i++) 

{  for  (j=0;j<=50;j++) 

{  u[n][i][j]  =  0.0; 

v[n][i][j]  =  0.0; 
U[n][i][j]  =  0.0; 
V[n][i][j]  =  0.0; 

} 

} 

} 


/*  Loops  */ 

time=200; 
ind  =  0; 
n=l; 

/*  Boundary  Conditions  */ 
for  (i=l;i<=49;i++) 

{  j=i; 

y=1.0; 

f=exp(-(i-25)*(i-25)); 

U[n][i][j-l]=U[n][i][j+l]+V[n][i+l][j]-V[n][i-l]|j]; 

u[n][i][j-l]=u[n][i][j+l]-v[n][i-l]lj]+v[n][i+l]0]; 

II=-(2*h*f*y)+((P-2*Mu)*(u[n][i+l][j]-u[n][i-l][j])+P*v[n][i]U+l] 
+Q*(V[n][i][j+l]-U[n][i-l][j]+U[n][i+l][j])); 
JJ=(beta*2*h*f*y)+Q*(v[n]  [i]  [j+ 1  ]-u[n]  [i- 1  ]  [j]+u[n]  [i+ 1  ][j])+ 

R*(  V[n]  [i]  [j+ 1  ]-U[n]  [i- 1  ]  [j  ]+U[n]  [i+ 1  ]  [j]) ; 

v[n][i]0-l]=(H*R  -  Q*JJ)/(P*R  -  Q*Q); 

V[n][i][j-1]=(Q*II  -  P*JJ)/(Q*Q  -  P*R); 

} 

while  (ind  <  time) 

{ 

for  (i=l;i<=49;i++) 

{  for  (j=l;j<=49;j++) 

{ 

Aa=u[n][i+l][j]; 
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Bb=u[n][i][j]; 

Cc=u[n][i-l][j]; 

Dd=u[n-l][i][j]; 

Ee=u[n][i][j+1]; 

Ff=u[n][i][j-1]; 

Gg=v[n][i-l][j+l]; 

Hh=v[n][i-l][j-l]; 

Ii=v[n][i+l][j+l]; 

Jj=v[n][i+1]0-1]; 

Kk=U[n][i+l][j]; 

Ll=U[n]  [i]  [j] ; 

Mm=U[n][i-l][j]; 

Nn=U[n-l][i][j]; 

Oo=V[n][i+l][j+l]; 

Pp=V[n][i-l]U+l]; 

Qq=V[n][i-l][j-l]; 

Rr= V  [n]  [i+ 1  ]  U  - 1  ] ; 

aA=u[n][i+l][j+l]; 

bB=u[n][i-l][j+l]; 

cC=u[n][i-l][j-l]; 

dD=u  [n]  [i+ 1  ]  [j  - 1  ] ; 

eE=v[n][i][j+l]; 

fF=v[n][i][j]; 

gG=v[n][i][j-l]; 

hH=v[n- 1  ]  [i]  [j] ; 

il=v[n][i+l][j]; 

jj=v[n][i-l][j]; 

kK=U  [n][i+l][j+l]; 

lL=U[n][i-l][j+l]; 

mM=U[n][i-l]|j-l]; 

nN=U[n][i+l][j-l]; 

oO=V[n][i][j+l]; 

pP=V[n][i][j]; 

qQ=V[n]  [i]  [j- 1  ] ; 

rR=V[n-l][i][j]; 


EE=FEE(Aa,Bb,Cc,Dd,Ee,Ff,Gg,Hh,Ii,Jj,Kk,Ll,Mm,Nn,Oo,Pp,Qq,Rr); 

FF=FFF(Aa,Bb,Cc,Dd,Ee,Ff,Gg,Hh,Ii,Jj,Kk,Ll,Mm,Nn,Oo,Pp,Qq,Rr); 

GG=FGG(aA,bB,cC,dD,eE,fF,gG,hH,iI,jJ,kK,lL,mM,nN,oO,pP,qQ,rR); 

HH=FHH(aA,bB,cC,dD,eE,fF,gG,hH,iI,jJ,kK,lL,mM,nN,oO,pP,qQ,rR); 

/*  Calculating  the  n+1  values  */ 
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u  [n+ 1  ]  [i]  [j]=(BB  *EE- A  A*FF)/(BB  *CC+AA*DD) ; 
v[n+l][i][j]=(BB*GG-AA*HH)/(BB*CC+AA*DD); 
U[n+l][i][j]=(DD*EE+CC*FF)/(BB*CC+AA*DD); 
V[n+l][i][j]=(DD*GG+CC*HH)/(BB*CC+AA*DD); 

} 

} 

for  (i=0;i<=50;i++) 

{  for  (j=0;j<=50;j++) 

{ 

/*  Updating  the  grid  for  next  time  level  */ 

u[n- 1]  [i]  [j]=u[n]  [i]  [j] ; 
u[n][i][j]=u[n+l][i][j]; 
v[n- 1  ]  [i]  [j]=v[n]  [i]  [j] ; 
v[n][i][j]=V[n+l][i][j]; 

U[n-l][i][j]=U[n][i][j]; 

U  [n]  [i]  [j]=U[n+ 1  ]  [i]  [j  ] ; 

V[n-l][i][j]=V[n][i][j]; 

V[n][i][j]=V[n+l][i][j]; 

} 

} 

for  (i=l;i<=49;i++) 

{ 

sol[ind]  [i]=sqrt(v[2]  [i]  [  1  ]  *  v[2]  [i]  [  1  ]+u[2]  [i]  [  1  ]  *u[2]  [i]  [  1  ]); 
flu[ind]  [i]=sqrt(  V[2]  [i]  [  1  ]  *  V[2]  [i]  [  1  ]+U  [2]  [i]  [  1  ]  *U[2]  [i]  [  1  ]) ; 

} 


ind  =  ind  +  1 ; 

} 

for  (tt=0;tt<time;tt++) 

{ 

for  (i=0;i<=50;i++) 

printf("\n%d\t%d\t%g\t%g",tt,i,sol[tt][i],flu[tt][i]); 

} 

/*for  (i=l;i<=49;i++) 

{ 
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for  (j=l;j<=49;j++) 

{ 

printf("\n%d\t%d\t%g\t%g\t%g\t%g",i,j  ,u[n]  [i]  [j]  ,v[n]  [1]  [j] , 
U[n][i][j],V[n][i][j]); 

} 

}*/ 

return  0; 

} 


/*  Declared  Function  Definition  */ 

float  FEE(float  Aa,float  Bb,float  Cc, float  Dd,  float  Ee, float  Ff, float  Gg, 
float  Hh, float  Ii,float  Jj, float  Kk, float  LI, float  Mm, float  Nn, 
float  Oo, float  Pp, float  Qq, float  Rr) 

{ 

a  =  ( 1  -beta+(beta*Kr/Kf)-(Kb/Kr)) ; 

Q  =  (beta*(l-beta)*Kr-beta*Kb)/a; 

R  =  beta*beta*Kr/a; 

P  =  (( 1  -beta)*((  1  -beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a  +  4*  Mu/3; 

b  =  N*beta*beta/K; 

rhol2  =  rhof*beta*(l -alpha); 

rholl  =  (l-beta)*rhor-rhol2; 

rho22  =  rhof*beta-rhol2; 

ans  =  Q/Mu*((Aa-2*Bb+Cc)/(h*h)+(Ii-Gg+Hh-Jj)/(4*h*h))  + 

R/Mu*((Kk-2*Ll+Mm)/(h*h)+(Oo-Pp+Qq-Rr)/(4*h*h))  - 
b*Fw*w*((Dd-Nn)/(2*delT*Mu))  - 

rho22*(P-Mu)/(Mu*rhol2)*((Aa-2*Bb+Cc)/(h*h)+(Ii-Gg+Hh-Jj)/(4*h*h))  - 

rho22/rho  1 2  *  ( Aa-4  *Bb+Cc+Ee+Ff)/(h*h)  - 

rho22*Q/(rhol2*Mu)*((Kk-2*Ll+Mm)/(h*h)+(Oo-Pp+Qq-Rr)/(4*h*h))  - 
rho22*b*Fw*w*(Dd-Nn)/(rhol2*2*delT*Mu)  - 
(Dd-2*Bb)*w*w*(rhol2-((rho22*rholl)/rhol2))/(delT*delT*Mu); 

return  ans; 

1 

float  FFF(float  Aa, float  Bb,float  Cc, float  Dd,float  Ee, float  Ff,float  Gg, 
float  Hh,float  Ii,float  Jj, float  Kk,float  LI, float  Mm, float  Nn, 
float  Oo, float  Pp,float  Qq,float  Rr) 

{ 

a  =  ( 1  -beta+(beta*Kr/Kf)-(Kb/Kr)); 

Q  =  (beta*  ( 1  -beta)*  Kr-beta*  Kb)/a; 

R  =  beta*beta*Kr/a; 


38 


P  =  (( 1  -beta)*((  1  -beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a  +  4*Mu/3; 

b  =  N*beta*beta/K; 

rhol2  =  rhof*beta*(  1-alpha); 

rholl  =  (l-beta)*rhor-rhol2; 

rho22  =  rhof*beta-rhol2; 

ans  =  Q/Mu*((Aa-2*Bb+Cc)/(h*h)+(Ii-Gg+Hh-Jj)/(4*h*h))  + 

R/Mu*((K3t-2*Ll+Mm)/(h*h)+(Oo-Pp+Qq-Rr)/(4*h*h))  - 
b*Fw*w*((Dd-Nn)/(2*delT*Mu))  - 

rhol2*(P-Mu)/(rholl*Mu)*((Aa-2*Bb+Cc)/(h*h)+(Ii-Gg+Hh-Jj)/(4*h*h))  - 
rho  1 2/rho  1 1  *(  Aa-4*Bb+Cc+Ee+Ff)/(h*h)  - 

rhol2*Q/(rholl*Mu)*((Kk-2*Ll+Mm)/(h*h)+(Oo-Pp+Qq-Rr)/(4*h:,:h))  - 
rhol2*b*Fw*w*(Dd-Nn)/(rholl*2*delT*Mu)  - 
(Nn-2*Ll)*w*w*(rho22-((rhol2*rhol2)/rhol  l))/(delT*delT*Mu); 
return  ans; 

} 

float  FGG(float  aA, float  bB, float  cC, float  dD,float  eE, float  fF, float  gG, 
float  hH,float  iI,float  j  J,float  kK,float  lL,float  mM, float  nN, 
float  oO, float  pP, float  qQ,float  rR) 

{ 

a  =  ( l-beta+(beta*Kr/Kf)-(Kb/Kr)); 

Q  =  (beta*(l-beta)*Kr-beta*Kb)/a; 

R  =  beta*beta*Kr/a; 

P  =  (( 1  -beta)*((  1  -beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a  +  4*Mu/3; 

b  =  N*beta*beta/K; 

rho  12  =  rhof*beta*(l -alpha); 

rholl  =  ( 1  -beta)*rhor-rho  12; 

rho22  =  rhof*beta-rhol2; 

ans  =  Q/Mu*((aA-bB+cC-dD)/(4*h*h)+(eE-2*fF+gG)/(h*h))  + 

R/Mu*((kK-lL+mM-nN)/(4*h*h)+(oO-2*pP+qQ)/(h*h))  - 
b*Fw*w*(hH-rR)/(2*delT*Mu)  - 

rho22*(P-Mu)/(rhol2*Mu)*((aA-bB+cC-dD)/(4*h*h)+(eE-2*fF+gG)/(h*h))  - 
rho22/rho  1 2*(iI-4*fF+j J+eE+gG)/(h*h)  - 

rho22*Q/(rhol2*Mu)*((kK-lL+mM-nN)/(4*h*h)+(oO-2*pP+qQ)/(h*h))  - 
rho22*b*Fw*w*(hH-rR)/(rho  1 2*2*delT*Mu)  - 
(hH-2*fF)*w*w*(rhol2-((rho22*rholl)/rhol2))/(delT*delT*Mu); 
return  ans; 

} 


float  FHH(float  aA,float  bB,float  cC,float  dD,float  eE,float  fF,float  gG, 
float  hH, float  il, float  j  J,float  kK, float  1L, float  mM,float  nN, 
float  oO,float  pP, float  qQ, float  rR) 

{ 
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a  =  ( 1  -beta+(beta*Kr/Kf)-(Kb/Kr)); 

Q  =  (beta*  (1 -beta)  *Kr-beta*Kb)/a; 

R  =  beta*beta*Kr/a; 

P  =  (( 1  -beta)*((  1  -beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a  +  4*Mu/3; 

b  =  N*beta*beta/K; 

rhol2  =  rhof*beta*(  1-alpha); 

rhol  1  =  (l-beta)*rhor-rhol2; 

rho22  =  rhof*beta-rhol2; 

ans  =  Q/Mu*((aA-bB+cC-dD)/(4*h*h)+(eE-2*fF+gG)/(h*h))  + 

R/Mu*((kK-lL+mM-nN)/(4*h*h)+(oO-2*pP+qQ)/(h*h))  - 
b*Fw*w*(hH-rR)/(2*delT*Mu)  - 

rho  1 2*(P-Mu)/(rho  1 1  *Mu)*((aA-bB+cC-dD)/(4*h*h)+(eE-2*fF+gG)/(h*h))  - 
rhol2/rhol  l*(iI-4*fF+jJ+eE+gG)/(h*h)  - 

rhol2*Q/(rhol  l*Mu)*((kK-lL+mM-nN)/(4*h*h)+(oO-2*pP+qQ)/(h*h))  - 
rho  1 2*b*Fw*w*(hH-rR)/(rho  1 1  *2*delT*Mu)  - 
(rR-2*pP)*w*w*(rho22-((rhol2*rhol2)/rhol  l))/(delT*delT*Mu); 
return  ans; 

} 
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